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ABSTRACT 

We employ high-resolution iV-body simulations of isolated spiral galaxy models, from low-amplitude, 
multi-armed galaxies to Milky Way-like disks, to estimate the vertical action of ensembles of stars 
in an axisymmetrical potential. In the multi-armed galaxy the low-amplitude arms represent tiny 
perturbations of the potential, hence the vertical action for a set of stars is conserved, although after 
several orbital periods of revolution the conservation degrades significantly. For a Milky Way-like 
galaxy with vigorous spiral activity and the formation of a bar, our results show that the potential 
is far from steady, implying that the action is not a constant of motion. Furthermore, because of the 
presence of high-amplitude arms and the bar, considerable in-plane and vertical heating occurs that 
forces stars to deviate from near-circular orbits, reducing the degree at which the actions are conserved 
for individual stars, in agreement with previous results, but also for ensembles of stars. If confirmed, 
this result has several implications, including the assertion that the thick disk of our Galaxy forms 
by radial migration of stars, under the assumption of the conservation of the action describing the 
vertical motion of stars. 

Subject headings: galaxies: kinematics and dynamics - Galaxy: disk - Galaxy: evolution - stars: 
kinematics and dynamics 


1. INTRODUCTION 

A stellar system can be fully modeled by a phase-space 
density distribution function (DF), whose evolution is de¬ 
scribed by solutions of the Boltzmann equation (Binney 
& Tremaine 2008). Is also a well-known result that for 
steady-state potentials, the DF describing the galaxy can 
be represented as a function of the integrals of motion 
of the system (Jeans 1919). There are several choices 
one can make to use then as coordinates of the DF, e.g. 
energy and angular momentum (Eddington 1915). In¬ 
tegrated orbits in axisymmetric potentials indicate the 
existence of a non-classical isolating integral of motion 
(Ollongren 1962), which is closely related to the actions 
of the system (Binney & Spergel 1984). 

Besides being a natural choice to describe the evolution 
in the phase-space of collisionless stellar systems, the ac¬ 
tions are also adiabatic invariants. This is a very appeal¬ 
ing property of spiral galaxies, where the high frequency 
of the vertical motion, compared to in-plane evolution, 
makes the action describing the z motion a prime candi¬ 
date for being a conserved quantity, even in the presence 
of spiral arms. This assumption has been tested by Sol¬ 
way et al. (2012), who showed that for isolated galaxies 
the vertical action is not a constant of motion of individ¬ 
ual stars, and is only conserved in average for samples of 
stars, with an intrinsic dispersion of ^ 20%. 

It is known that recurring spiral arm-activity in a 
galactic disk causes stars to diffuse radially over time in 
a manner that largely preserves the overall structure of 
the disk. Sellwood & Binney (2002) used bi-dimensional 
stellar disks to show that stars near the corotation radius 
of a transient spiral pattern change angular momentum 
and move to new radii inwards or outwards with negligi¬ 
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ble increase of random motion, causing a process termed 
radial migration. Following studies confirmed that this 
process occurs in 3D simulations of galactic disks either 
in isolation or in a cosmological context (Roskar et al. 
2008; Minchev et al. 2011; Grand et al. 2012; Loebman 
et al. 2011; Kubryk et al. 2014; Halle et al. 2015; Grand 
et al. 2015). 

If radial mixing in galactic disks occurs, this process 
has interesting implications for the spread in the age- 
metallicity relation observed in the solar neighborhood 
(Edvardsson et al. 1993). If stellar migration never oc¬ 
curred, the stars in the solar vicinity would have metallic- 
ities ranging from zero to the present-day value (Serenelli 
et al. 2009) with a significant correlation between the age 
of the stars and their metallicities. However this correla¬ 
tion has never been found in the current data, raising the 
possibility that the stars migrated from their birth radii 
(Wielen et al. 1996). Because the radii of stars at the 
birth depends on their specific angular momentum, stel¬ 
lar migration occurs if there is a change of the angular 
momentum of the stars. Changes in the angular mo¬ 
mentum for individual stars can be induced by scatter¬ 
ing with giant molecular clouds (Spitzer & Schwarzschild 
1953; Mihalas & Binney 1981) or transient spiral struc¬ 
tures (Jenkins & Binney 1990; Jenkins 1992), interfer¬ 
ence of coherent long-lived modes and density waves with 
a bar (Roskar et al. 2012), or by interactions with satel¬ 
lite galaxies (Kazantzidis et al. 2008; Villalobos & Helmi 
2008; Bird et al. 2012). 

In this paper we aim to study in detail the effect of 
migration on isolated disks, by using TV-body simula¬ 
tions for disk galaxies evolving in isolation, to achieve 
high-resolution; with the aim of better determining the 
contribution of the spiral waves to the heating of the stel¬ 
lar disk, to the radial migration of stars in the disk, and 
ultimately to the formation of a thick disk. 

This manuscript is laid out as follows. In Section 2 we 
introduce the TV-body simulations used in our work and 
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describe the method to measure the vertical action. In 
Section 3 we discuss to what extent this quantity is con¬ 
served during the evolution of the disk, and in Section 4 
we show how this can be used to study the formation of a 
thick disk. Finally in Section 5 we draw our conclusions. 

2. NUMERICAL PRELIMINARIES 
2.1. Simulations 

For this study we employed TV-body simulations of 
two galaxy models labeled as disk 1 and disk 2. Each 
galaxy consisted of a static dark matter halo and a live 
rotationally supported stellar disk of 5 x 10 6 particles. 
Because the spiral morphology in stellar disks depends 
in detail on the mass distribution of the galaxy and 
on the self-gravity of the disk, we selected the struc¬ 
tural parameters of the galaxies such that one galaxy 
developed a low-amplitude multi-armed disk (disk 1), 
whereas the other developed a spiral morphology more 
similar to a Milky Way (MW)-like galaxy (disk 2). This 
was achieved by choosing models with very different disk 
mass fractions within 2.2 scale-lengths: Moi s k/M tota \ = 
0.271 for disk 1 and 0.478 for disk 2, respectively. By 
varying the mass distribution of the disk within 2.2 scale- 
lengths, these models were designed to change the critical 
length scale parameter A cr i t = 47r 2 GE/ft 2 , where E is the 
disk surface mass density and k the epicycle frequency, 
and lead to a different spiral morphology (Toomre 1981; 
Sellwood & Carlberg 1984; Carlberg & Freedman 1985; 
D’Onghia 2015). 

The parameters describing each galaxy component 
were assumed to be independent and both models were 
constructed similarly to the approach described in (Hern- 
quist 1993; Springel et al. 2005). The dark matter mass 
distribution was modeled assuming the Hernquist profile 
(Hernquist 1990): 


Phalo(r) 27ri?3 alo (r/i? halo )(l + r/i? halo )3 ’ (1) 

with Mh a io being the total halo mass and Rhaio the scale 
radius. The effect of the dark halo is modeled assuming a 
static potential, this is of course a limitation of our model 
since the system never completely relaxes, however, the 
complication of adding an alive component comes with 
the price of including further scattering to the particles 
of the disk. The stellar disk is modeled in the initial 
conditions as a thin exponential surface density profile 
of scale-length i?disk and total mass M^isk- The vertical 
mass distribution of the stars in the disk is specified by 
the profile of an isothermal sheet with a radially constant 
scale height Zdisk- 


Pdisk(-^5 


serh 2 _JL 
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Model disk 2 included a static potential for the bulge 
described by the Hernquist model with mass Mb u i ge and 
scale-length Rbuige- We aim to build this galaxy model 
with a flat rotation curve consistent with the observa¬ 
tions of the MW within two optical radii (R < 14 kpc) 
(Reid et al. 2014; Bovy et al. 2009). We noted that this 
requirement is satisfied by adopting a Hernquist profile 
for the dark halo with a large scale radius, and a large 


TABLE 1 

Structural properties of the galaxy models 



disk 1 

disk 2 

M disk (10 1U M Q ) 

1.905 

4.000 

#disk (kpc) 

3.130 

2.500 

^disk/ A-disk 

0.100 

0.150 

M halo (10 12 M 0 ) 

0.933 

9.500 

#halo (kpc) 

29.775 

130.0 

-^Viisk/ALotal G'^-^disk) 

0.271 

0.479 

M bulge (10 10 M 0 ) 


1.400 

-^bulge (kpc) 


0.350 


total mass. However, the dynamical range of interest of 
this study is a few optical radii, thus the mass distri¬ 
bution outside this radius does not affect our analysis. 
Table 1 summarizes the structural parameters adopted 
for each galaxy model in this study. 

Fig. 1 shows the azimuthally averaged radial profiles 
for the total circular velocity Wire (stars and dark mat¬ 
ter), the Toomre parameter (^Toomre of the stellar disk, 
the disk surface density E and its vertical velocity dis¬ 
persion (j z . The profiles are displayed at four different 
times, spanning 6 Gyr of evolution. In what follows we 
use the notation t n to refer to the time of integration 
in our simulations in Gyr, with to thus labeling the the 
initial time, and to corresponding to 6 Gyr. 

We note that in the multi-armed galaxy (disk 1 
model), QToomre increases with time, as a consequence 
of the increase of the radial velocity dispersion. This is 
expected, since in-plane heating occurs as a consequence 
of the spiral structure activity, while the vertical struc¬ 
ture remains unchanged. 

This is, however, not the case for our disk 2 run. A 
bar instability grows at the center of the galaxy after 
~ 3 Gyr. Despite the flat rotation curve, the surface 
mass density of the disk changes considerably with time 
and the vertical velocity dispersion increases by ^ 15% 
in the shown radial range. 

In Fig. 2 a sequence in time of the face-on views of the 
stellar disk of both models of galaxies are displayed. The 
colored bar indicates the values for the surface density of 
the disk. As expected in the context of swing amplifi¬ 
cation theory, the two models show very different spiral 
morphologies. Our disk 1 model grows a large number 
of low-amplitude spiral arms that self-perpetuate for at 
least 6 Gyr (D’Onghia et al. 2013). After 1 Gyr disk 2 
already developed four high-amplitude spiral arms at two 
scale-lengths (^ 5 kpc), in remarkable agreement with 
the number of arms predicted for a disk with 50% disk 
fraction (D’Onghia 2015; Pettitt et al. 2015). Later in 
time this disk model develops a bar (Merritt & Sellwood 
1994; Martinez-Valpuesta et al. 2006; Athanassoula & 
Martinez-Valpuesta 2008; Yurin & Springel 2015), which 
becomes the dominant non-axisymmetric feature after 
rsj 4 Gyr of evolution. 

Note that the efficiency of the stellar radial migration 
by recurrent spiral activity in the multi-armed galaxy 
(disk 1 model) has been recently studied by Vera-Ciro 
et al. (2014). As first shown in Sellwood & Binney (2002), 
the authors confirmed that spiral activity in a low-mass 
disc causes stars to diffuse radially, with stars near coro¬ 
tation of a spiral pattern exchanging angular momentum 
and moving to a new radius without adding random mo¬ 
tion (see also Daniel & Wyse 2015). Indeed, the mech- 
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Fig. 1. — Radial profiles of total circular velocity, Fb r , Toomre 
parameter, QToomre? surface mass density, S, and vertical velocity 
dispersion, cr z , for the two studied galaxy models: disk 1 (left) and 
disk 2 (right). t n refers to the time of integration in our simula¬ 
tions in Gyr, with to being the the initial time and tQ corresponding 
to 6 Gyr. 

anism of scattering of stars at corotation predicts a lack 
of significant disk heating and leaves the surface density 
profile unchanged. These two features are reproduced in 
our disk 1 model as shown in Fig. 2 (left panels). In 
the MW-like galaxy model (disk 2) the bar formation 
causes greater angular momentum changes than those 
of a multi-armed spiral structure (Minchev & Famaey 
2010). Furthermore, because the bar persists after its 
formation, the associated angular momentum changes 
might differ from those caused by recursive spiral arms. 
Therefore the surface density profile of the disk is ex¬ 
pected to change. Some disk heating is also expected to 
occur with time. 

2.2. Measuring the vertical action 

The action-angle variables can be used to describe the 
evolution of orbits in static potentials. The actions, j, 
are constants of motion in the unperturbed field, or when 
the potential changes slowly so that they are adiabatic 
invariants (Arnold 1978; Goldstein et al. 2002). The ver¬ 
tical action is defined as: 

i* = i l dz dv * = 4 / v '*> ( 3 ) 


Fig. 2. — Time sequence of density projections of the face-on 
views of stellar disks presented in this study. The disk 1 simula¬ 
tion develops a large number of low-amplitude spiral structure that 
perpetuate for at least 6 Gyr (left panels); disk 2 develops approx¬ 
imately a four-fold rotational symmetry and becomes bar-unstable 
after ~ 3 Gyr. 

nearly axisymmetric potential. Solway et al. (2012), for 
instance, calculated the quadrature in Eq. (3) by mea¬ 
suring the area enclosed by the orbit in the z — v z plane. 
We adopted here the prescription as described in Ap¬ 
pendix A. This method allows us to identify and reject 
stars trapped in resonances, a case where the adiabatic 
invariance breaks (Pfenniger 1984; Pfenniger & Friedli 
1991). 

In Fig. 3 we show an example for a random particle in 
our disk 1 simulation at two different times. The green 
points show the surface of section of the integrated orbit 
and the red line is the connected loop used to calculate 
Eq. (3). The blue line represents the zero-velocity curve, 
and gray points are other orbits, shown here to depict 
the structure of the phase-space in our runs. The area 
defined by the red-loop in Fig. 3 is the vertical action j z , 
for this particular example, the vertical action measured 
at t\ is a factor of ^ 80 larger than it is at time £3, this 
is clear evidence that the value of j z is not conserved for 
individual orbits for slow-varying, low-amplitude pertur¬ 
bations, such as the ones developed in our disk 1 run. 


where z and v z describe the position and velocity along 
the orbit, respectively. There are a number of ways to 
estimate of the vertical action of a star particle in a 


3. A CONSERVED QUANTITY? 

Our previous calculations demonstrated that the ver¬ 
tical action j z is not a conserved quantity for individ- 
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Fig. 3. — Surface of section for a star particle measured at two 
different times. In the left panel the simulation has already devel¬ 
oped spiral arms, the right panel is the surface of section after 2 
Gyr. In this specific example, R gu i(ti) ~ R gu i(t3) indicates that 
the particle has not migrated significantly. However, its vertical 
action has changed by a factor ~ 80 ( Sj z = 4.4). For each case the 
solid red line displays the closed path used to define the action and 
the blue line indicates the zero-velocity curve. 
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Fig. 4. — Probability distribution of fractional changes of the 
vertical action j z , as a function of changes in guiding radius R gu i 
for particles selected at 1 Gyr with R gu i = 5.00 d= 0.25 kpc. The 
distribution of fractional changes in j z does not depend on how 
efficiently particles migrate. The median and lcr equivalent disper¬ 
sion are represented with solid lines in the left panels. The right 
panels show the marginalized distribution of changes in the vertical 
action. 


ual stars. This result agrees with findings of previous 
studies (Solway et al. 2012). We investigated this prop¬ 
erty in our disk 1 and disk 2 simulations. After 1 Gyr 
of evolution an ensemble of stars with guiding radius: 
R g ui = 5.00 d= 0.25 kpc was selected and the relative 
changes in the vertical action were computed at a later 
time (See Appendix A). This radius is close to two disk 
scale-lengths in both galaxy models. The vertical action 


changed significantly over a period of time as shown in 
Fig. 3. Thus, its change from time t to a later time t' > t 
was quantified in the following way: 

Sj z = ln^%. (4) 

3z(t) K 1 

In addition, we noted that over the same period of time 
a significant radial migration occurs, with a consequent 
change of the guiding radii. Thus, a similar definition 
was also applied to SR gu p Fig. 4 (top panel) displays the 
outcome of our study for the multi-armed galaxy disk 

1 over a period of 2 Gyr. In agreement with previous 
studies, our findings suggest the vertical action is a con¬ 
served quantity on average for an ensemble of stars. Note 
that the modest fractional changes of 5j z in our simula¬ 
tion match the values reported by Solway et al. (2012) 3 
Furthermore, Fig. 4 indicates that any variation in the 
vertical action, 5j z , is independent of the radial migra¬ 
tion, measured by changes in guiding radii, SR gu [. Next, 
we computed the median of the fractional change in the 
vertical action n(5j z ) and the lcr equivalent dispersion 
around the median a(Sj z ). We obtained n(bj z ) = 0, 
with standard deviation cr(5j z ) = 0.28 (displayed as the 
solid black line in Fig. 4), also in agreement with the 
values reported in Solway et al. (2012). 

The same analysis performed on the more massive disk 

2 produces different results as shown in Fig. 4 (bottom 

panel). The median of fractional changes in the vertical 
action is = 0.09 with a dispersion a(Sj z ) = 0.64, 

significantly larger than the case of the multi-armed 
galaxy. Thus, changes of angular momentum occurred 
with the formation of the bar, with significant disk heat¬ 
ing, combined with the strong variation of the potential 
of the disk, do not conserve the vertical action at the 
same level of accuracy. 

To demonstrate this unexpected outcome, the evolu¬ 
tion by time of the median fjb(5j z ) and its dispersion 
cr(5j z ) of fractional changes in the vertical action are dis¬ 
played in Fig. 5 for the two galaxy models. Ensembles of 
stars are selected at after 1 Gyr of evolution (t = t±) with 
guiding radii R gu i = 3, 5, 8 kpc and followed over the evo¬ 
lution time of the disk. While in the multi-armed disk, 
even after 40 orbital periods the average variation in the 
vertical action is never larger than 1%, in the MW-like 
disk galaxy the median reaches values of 15-20% after 20 
orbital periods. 

Note that in Fig. 5 the dispersion around the median 
steadily increases as a function of time for both mod¬ 
els (left panels), suggesting that the vertical action of an 
ensemble of stars tends to be less conserved as time pro¬ 
gresses. This results holds for both the multi-armed disk 
and the MW-like disk. 

We also checked whether the vertical heating of the 
disk, measured by the vertical velocity dispersion of the 
ensemble of stars, could produce changes in the average 
vertical action. Fig. 5 plots the median n(5j z ) and dis¬ 
persion cr(Sj z ) of fractional changes in the vertical action 
against the corresponding variation in the velocity dis¬ 
persion of the vertical component. 

To test for the origin of the changes in the average ac¬ 
tion we note that the vertical energy of samples of stars 

3 \nx w x — 1, for x ~ 1. 
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Fig. 5. — Left panels: evolution of average changes in the vertical 
action (top panels) and its la equivalent dispersion (bottom pan¬ 
els) as a function of time. Right panels: the same quantities as in 
the left panels are plotted against the changes in vertical velocity 
dispersion for ensembles of stars selected at different radii. 

proportionally depends on the vertical velocity dispersion 
(j z ) rsj (j\ (cf. Eq. (7)), and so if the average vertical ac¬ 
tion increases as a consequence of the vertical heating in 
the disk, then 5j z = 2 5a z . The multiplying factor in this 
approximation is the vertical epicycle frequency, which 
remains mostly constant in time (See Fig. 1). Fig. 5 
shows that this is not the case, and that the increase 
in dispersion of the changes in the vertical action is not 
merely determined by the vertical heating of the disk. 

4. IMPLICATIONS FOR THE VERTICAL STRUCTURE 

What does the conservation of the vertical action im¬ 
ply for the structure and evolution of the stellar disk? It 
has been argued that under the assumption of the conser¬ 
vation of the vertical action j z stars migrating outwards 
in the disk tend to move far from the mid-plane and 
populate the thick disk (Roskar et al. 2008; Schonrich & 
Binney 2009). We address this question here. 

The epicycle approximation regulates the motion of 
stars on near-circular orbits. In the approximation the 
vertical motion of a star is also decoupled from the in¬ 
plane motion. The motion perpendicular to the plane is 
then described by the Hamiltonian (Binney & Tremaine 
2008): 

Hz = \v 2 z + X -v 2 z 2 . (5) 

For an ensemble of star particles at radius i?, the av¬ 
erage vertical energy is 

(H z )= 1 -a 2 z + 1 -u 2 (z 2 ), (6) 

where cr z is the vertical velocity dispersion of the ensem¬ 
ble. Note that o z in Eq. (6) does not necessarily refer 


Fig. 6 . — Top: average vertical action as a function of radius, the 
solid lines are the values calculated using the method described in 
the previous section (Eq. (3)) whereas the dashed lines correspond 
to the epicycle approximation (Eq. (6)). The difference between 
these two values is ~ 15%. Bottom: thickness of the disk calculated 
as the RMS values of the ^-coordinates (solid), and by solving 
Eq. (7). 

to the value of velocity dispersion of all stars with home 
radii at R. Indeed, Eq. (6) holds also for any popula¬ 
tion of stars whose guiding radii have changed over time 
(e.g. stars migrating from inner or outer locations within 
the disk) and that are used to calculate it. Thus, in the 
epicycle approximation, the vertical action is the follow¬ 
ing j z = H z /v (Binney & Tremaine 2008): 

{jz,epi) = 2F = T(o-2 + jy 2 a 2 h 2 ), (7) 

where, a 2 h 2 = (z 2 ), ( z 2 ) is the second moment of the 
vertical component, and a = tt/\/12. With this nor¬ 
malization, the parameter h is equivalent to the disk 
scale height Zdisk, h = ^disk, when the second moment 
is calculated for all stars sampled from the sech 2 z/zdi s k 
model used for the vertical density of our simulations (cf. 
Eq. (2)). 

Once the average vertical action and vertical velocity 
dispersion are computed for an ensemble of stars, Eq. (7) 
is solved for the thickness h of that population of stars. 

Fig. 6 displays the average value of the vertical action 
j z as a function of the galactic radius at three different 
times for both galaxy models (top panels). The solid 
lines show the average vertical action calculated for a set 
of stars using the algorithm described in Section 2. The 
dashed lines illustrate the predicted vertical action using 
the epicycle approximation for that set of stars. Clearly, 
the epicycle approximation predicts values of the vertical 
action 15% higher than the estimates obtained with the 
numerical procedure. A similar trend is observed in the 
numerical experiments described in Solway et al. (2012). 
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Fig. 7.— Scale height h profile for particles in circular orbits 
that migrate from R gu i(ti) = 5.00d= 0.25 kpc (left) and -Rgui(U) = 
8.00 d= 0.40 (right) as a function of their final guiding radii for the 
the multi-arm simulation (top) and the MW-like model (bottom). 
The scale height profile of the whole disk at 1 3 is shown by the solid 
blue line. In green we show the thickness profile of the migrators 
at t\ and in red the thickness at £ 3 . The orange dashed line is 
the result of using Eq. (7) and assuming that the vertical action is 
conserved (j Zje p i(ti)> = O’^epife))- 

Then, Eq. (7) is used to compute the thickness of the 
disk h at each radius, and the values obtained for the 
both galaxy models are shown in Fig. 6 (bottom pan¬ 
els). A similar discrepancy is reflected in the estimate of 
the scale height at different radii. The values derived in 
the epicycle approximation are 5% larger than the values 
obtained by evaluating Eq. (3) as previously described. 

In both simulations, after the disks evolved for 1 
Gyr (ti), we selected a set of stars with guiding radii 
.Rg U i = 5.00 zb 0.25 kpc (left panels in Fig. 7) and 
.Rg U i = 8.00 =b 0.40 kpc (right panels), we follow these 
stars for two billion years longer (£ 3 ) and calculated their 
scale height h profile as a function of their final position 
in the disk. This would show if a subsample of migrating 
stars can change their thickness h as they migrate. 

We only select stars that are in nearly circular orbits 
L z /L c i rc > 0.97, where L z is the vertical angular mo¬ 
mentum and L c - irc is the angular momentum that a star 
has, should it move in a circular orbit with the same 
energy. By adding this constraint at both times £1 and 
£3 we ensure that the sample contains stars that migrate 
while preserving their circular orbits (Sellwood & Binney 
2002 ). 

For this set of particles the average vertical action in 
the epicycle approximation is also j z is computed at 
both times t\ and £ 3 . The thickness at both times is 
also included, in general stars that move toward the cen¬ 
ter tend to reduce their scale height, whereas stars that 
migrate outwards increase it. For comparison we also 
show in blue the scale height of the disk as a function 
of radius, calculated with the normalized RMS value of 


the ^-coordinates of all particles in the disk. This value 
matches the one reported as 2 Aisk in table 1 (blue solid 
line). 

The bias reported in (Vera-Ciro et al. 2014) is observed 
for both disks. This explains the fact that for stars that 
migrate outwards, the slight increase in their vertical 
scale height is not enough to create a thicker component 
at their final location. We also included in Fig. 7 the 
predicted thickness h at final £3 under the assumption 
that the action is conserved. It can clearly be seen that 
for our disk 1 no thick distribution of stars is either ob¬ 
served or predicted even when the average vertical action 
is clearly conserved. 

In the MW-like galaxy, disk 2, a bar grows in the in¬ 
ner parts of the disk, in qualitative agreement with sim¬ 
ulations of Brunetti et al. (2011). As a consequence of 
the bar formation and the vigorous spiral activity there 
are not only significant changes of angular momentum, 
but also a considerable in-plane heating that forces stars 
in this simulation to deviate from a near-circular orbits. 
There are good reasons to believe that in this case the en¬ 
ergy of vertical motion is clearly not decoupled from that 
of the radial part of the motion. Indeed, the epicycle ap¬ 
proximation is a poor description of the motion of most 
particles, for which the radial and vertical oscillations are 
neither harmonic nor are the energies of the two oscilla¬ 
tions decoupled. Furthermore, the actions are constants 
of motion under the conditions of slow changes of the 
potential. In the multi-armed galaxy the low-amplitude 
arms represent tiny perturbations of the potential, hence 
the vertical action might still be conserved, although af¬ 
ter several orbital periods of revolution the vertical ac¬ 
tion is not conserved at the same level of accuracy. How¬ 
ever, in the MW-like galaxy the high-amplitude spiral 
structure and the bar formation seem to compromise the 
validity of these conditions. 

5. SUMMARY AND CONCLUSIONS 

We have presented quantitative estimates of the verti¬ 
cal action of ensembles of stars in the presence of spiral 
activity. The vertical action of star particles has been 
evaluated in isolated disks of galaxies with different spi¬ 
ral morphologies. The two galaxy models consist of a 
low-mass multi-armed disk (disk 1) and a MW-like disk 
that develops a bar after ^3.5 Gyr of evolution (disk 

2). 

For the multi-armed disk galaxy, disk 1 , we find that 
the vertical action for an ensemble of stars is approxi¬ 
mately a conserved quantity, with a dispersion of frac¬ 
tional changes around its median of 20 %, in agreement 
with the results reported in the literature. However, as 
time progresses the action is conserved with less accu¬ 
racy, as indicated by the increase of the standard devi¬ 
ation up to the value of 50% after 6 Gyr of evolution. 
This also holds for a subset of stars that radially mi¬ 
grates by resonant scattering at corotation. These re¬ 
sults have implications for the vertical structure of the 
disk when stars radially migrate to the outer part of the 
disk. In this model stars migrating outwards are a heav¬ 
ily biased subset of stars with preferentially low vertical 
velocity dispersion. Thus, extreme migrators outwards 
are characterized by having a have small amplitude of 
their vertical excursion, independent of the conservation 
of the vertical action. 
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In the MW-like galaxy model, disk 2, the bar for¬ 
mation causes large angular momentum changes, with 
significant in-plane and vertical heating. Thus, in the 
presence of the formation of the bar and combined with 
significant disk heating, our results show that the verti¬ 
cal action cannot be assumed a constant of motion. In¬ 
deed, in this case the fractional changes in the estimates 
of the action for a set of stars progressively increase by 
time, with a standard deviation that approaches 80% af¬ 
ter 6 Gyr of disk evolution. If confirmed, this result casts 
doubts on the possibility of predicting the vertical struc¬ 
ture of the disk for a set of extreme migrators outwards, 
by assuming that their vertical action is conserved. 
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APPENDIX 

A. MEASURING THE VERTICAL ACTION 

The vertical action j z is defined as: 



/ 


d z dv z 



(Al) 


where z and v z describe the position and velocity along the orbit, respectively. If the potential is nearly axisymmetric 
this quadrature can be calculated using the algorithm described in (Solway et al. 2012). We adopted here a similar 
prescription as explained below: 


Vera-Ciro & D’Onghia 
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Fig. Al. — Surface of section for a star particle measured at two different times. In the left panel the simulation has already developed 
spiral arms, the right panel is the surface of section after 2 Gyr. In this specific example, Rg U i(ti) ~ -R gu i(^3) indicates that the particle has 
not migrated significantly. However, the particle moves away from a resonant orbit, a situation that is easily identified with our scheme. 

1 . At a given time t the gravitational potential <f>, its radial gradient d$/dR, and vertical gradient dQ/dz are 
computed on a polar grid (R, 0 , z) centered at the minimum of the (fixed) halo potential. The grid dimension 
is Nr x No x N z = 128 x 128 x 128, with the cylindrical coordinates ranging as follows: e$ < R < 30 kpc 
(logarithmically spaced), —i r < 0 < tt and —6 kpc < z < 6 kpc, where = 50 pc is the softening length. Note 
that these ranges span six times the scale-length and twenty times the scale height. The potential is numerically 
calculated at each node of the grid using a tree algorithm that includes up to the quadrupole contribution. We 
do not impose a symmetry about the galactic mid-plane, in fact, modest asymmetries that can develop over time 
are naturally captured by the grid. 

2 . At each point (R, z) the potential is calculated by averaging over the Nq = 128 azimuthal points: 


-i iyi 6 

— (A2) 

Ne U 

Equivalent expressions are derived for d&/dR and d&/dz. This approximation relies on the fact that the potential 
is nearly axisymmetric, this however may not be the case after if a bar grows at the center of the galaxy. 

Once the potential is computed at t time, for each star particle the following quantities are evaluated: the vertical 
angular momentum L z , the guiding radius R gu [ as the solution to the equation: 


L z — -Rguikcirc 


(A^g U i), 


(A3) 


and the vertical action j z , with the procedure described below. 

First, the position and velocity vectors are transformed to cylindrical coordinates (R, z). This point is integrated in 
the fixed potential defined by the grid for a time £ max = 10 S T Z , where T z = 2 tt/ is is the period of vertical oscillations, 
and is is the vertical epicycle frequency (Binney & Tremaine 2008) 


is 


2 


<9 2 $(R, z) 
dz 2 


z=0 


(A4) 


The vertical frequency can be easily calculated using five-point stencil derivatives applied to the d^/dz grid. 

The integration is performed using a symplectic integrator of order 4 (Candy & Rozmus 1991). With this choice, 
typical deviations in the energy are always smaller than 10 -6 and show no secular growth over the integration time 
£ max . This procedure yields a set of discrete points along the orbit of the form z^VR^Vz^i. 

During the integration of each star, every time the condition Ri < R gu [ < Ri+i , vrj > 0 is satisfied, the dynamical 
states are changed from {R(t), z(t),VR(t),v z (t)} to {t(R), z(R),vr(R),v z (R)} which is numerically integrated in the 
range Ri < R < R gu i- With this change of coordinates, the symplectic structure clearly breaks, therefore we have to 
use a different (non-conservative) integrator. In this case we use the explicit embedded Runge-Kutta Prince-Dormand 
( 8 , 9) method (Dormand & Prince 1980) implemented in the Gsl 4 library. After the exact intersection with the 
R = Rgui plane is found, the integration is resumed from Ri+\ with our symplectic method. 

In Fig. Al we show an example for a random particle in our disk 1 simulation. The green points show the set 
obtained with the procedure described above. The blue line represents the zero-velocity curve, and gray points are 
other orbits, shown here to depict the structure of the phase-space in our runs. At a first glance, the set in the left 
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panel resembles a simple closed structure in the (z,v z )~ plane, however, a careful inspection reveals a more complex 
structure. Indeed, the surface of section of this particle is composed of five resonant islands. At later times, however, 
the orbit traces a single closed structure (right panel). This transition is non-adiabatic and as such we would like to 
devise a way to identify these cases. 

The task of calculating the vertical action for a given star is then reduced to measuring the area enclosed by the set 
of points in the (z, r^-plane found above as the star moves in the galaxy. To achieve this, the points in the surface of 
section were connected with the following procedure: 

1. In the (z,v z )~plane a random point po = (zo, v z,o) is selected and its closest neighbor pi is found. For simplicity, 
both coordinates (z,v z ) are rescaled in the interval (—1,1). 

2. The closest N points to pi, {qk}k=x that have not been already included in a previous set are then selected. For 
each point qk the angles between the segments popi and piqu are defined. 

3. The point with minimum angle weighted by the distance is found and defined as the the new p\. We noted that 
in our implementation the weights ~ \piqk| 3 perform reasonably well. 

4. Steps 2 and 3 are repeated until po is the next point to be included in the cycle. 

5. Step 1 is then repeated until all points have been included in our analysis. 

This algorithm then yields a set of unconnected loops in the surface of section defined by R = R gu [ and vr> 0. In 
the cases where only one loop is present the action is just the area of such cycle and can be calculated using the popular 
formula based on the Green theorem. Otherwise, the particle is trapped in a resonance and it is no longer considered 
in our results. This implementation uses TV = 16 and rejects any cycle that has fewer points than this threshold. In 
addition, orbits for which the distance between the last and first point is larger than five times the average distance 
between the other points in the cycle are rejected. This choice avoids cases in which the algorithm ends before closing 
a loop, a circumstance that occurs for stars whose orbit takes them close to the zero-velocity curve. If the orbit is 
composed of several loops (e.g. left panel in Fig. 3) and the procedure fails to calculate the area for one of those loops, 
then the orbit is rejected altogether. 

Fig. A1 shows the results of the implementation of our algorithm to the particles of our simulations, as illustrated 
by the red line of the inset of the left panel. In our numerical experiments the procedure successfully computes the 
vertical action for ^ 85% of the particles located at a distance of R > 2 kpc from the disk center, and ^ 60% of 
particles orbiting at smaller radii. 





